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Abstract 

We study a contact process (CP) with two species that interact in a symbiotic manner. 
In our model, each site of a lattice may be vacant or host individuals of species A and/or 
B; multiple occupancy by the same species is prohibited. Symbiosis is represented by a 
reduced death rate, fj, < 1, for individuals at sites with both species present. Otherwise, 
the dynamics is that of the basic CP, with creation (at vacant neighbor sites) at rate A 
and death of (isolated) individuals at a rate of unity. Mean-field theory and Monte Carlo 
simulation show that the critical creation rate, \ c (p>), is a decreasing function of fi, even 
though a single-species population must go extinct for A < A c (l), the critical point of the 
basic CP. Extensive simulations yield results for critical behavior that are compatible with 
the directed percolation (DP) universality class, but with unusually strong corrections to 
scaling. A field-theoretical argument supports the conclusion of DP critical behavior. 
We obtain similar results for a CP with creation at second-neighbor sites and enhanced 
survival at first neighbors, in the form of an annihilation rate that decreases with the 
number of occupied first neighbors. 

PACS numbers: 05.50.+q,05.70.Ln,05.70.Jk,02.50.Ey,87.23.Cc 
Keywords: 
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I. INTRODUCTION 



Absorbing-state phase transitions have attracted much interest in recent decades, 
as they appear in a wide variety of problems, such as population dynamics, hetero- 
geneous catalysis, interface growth, and epidemic spreading Interest 
in such transitions has been further stimulated by recent experimental realizations 
3, Is]. 

The absorbing-state universality class associated with directed percolation (DP) 
has proven to be particularly robust. DP-like behavior appears to be generic for 
absorbing-state transitions in models with short-range interactions and lacking a 
conserved density or symmetry beyond translational invariance 3, [lo| - By contrast, 
models possessing two absorbing states linked by particle-hole symmetry belong to 
the voter model universality class fll |. 

The contact process (CP) pjj is probably the best understood model exhibiting 
an absorbing-state phase transition; it has been known for many years to belong to 
the DP class. The CP can be interpreted as a stochastic birth-and-death process 
with a spatial structure. As a control parameter (the reproduction rate A) is varied, 
the system undergoes a phase transition between extinction and survival. In this 
context it is natural to seek a manner to include symbiotic interactions in the CP. 
In the present work, this is done by allowing two CPs (designated as species A and 
B) to inhabit the same lattice. The two species interact via a reduced death rate, 
/i, at sites occupied by individuals of both species. (Aside from this interaction, 
the two populations evolve independently.) We find, using mean-field theory and 
Monte Carlo simulation, that the symbiotic interaction favors survival of a mixed 
population, in that the critical reproduction rate A c decreases as we reduce fi. Note 
that for A(/i) < A < A(l), only mixed populations survive; in isolation, either species 
must go extinct. 

In addition to its interest as a simple model of symbiosis, the critical behavior of 
the two-species CP is intriguing in the context of nonequilibrium universality classes. 
By analogy with the (equilibrium) n-vector model, in which the critical exponents 
depend on the number of spin components n, one might imagine that the presence 
of two species would modify the critical behavior. Using extensive simulations, we 



find that the critical behavior is consistent with that of directed percolation (DP), 
although with surprisingly strong corrections to scaling. An argument based on 
field theory supports the conclusion of DP scaling. We note that our result agrees 
with that of Janssen, who studied general multi-species DP processes [l^. Similar 
conclusions apply to a related model, a CP with creation at second-neighbor sites 
and enhanced survival at first neighbors, in the form of an annihilation rate that 
decreases with the number of occupied first neighbors. (In this case the two species 
inhabit distinct sub lattices.) 

The balance of this paper is organized as follows. In the next section we define 
the models and analyze them using mean-field theory. In Sec. Ill we present our 
simulation results, and in Sec. IV we discuss a field-theoretic approach. Sec. V is 
devoted to discussion and conclusions. 



II. MODELS AND MEAN-FIELD THEORY 

To begin we review the definition of the basic contact process. Following the 
usual nomenclature, we refer to an active site as being occupied by a "particle" and 
an inactive one as "vacant" . The CP 12| is a stochastic interacting particle system 



defined on a lattice, with each site i either occupied by a particle [&i(t) = 1], or 
vacant [(Ti(t) = 0]. Transitions from Oi = 1 to Oi = occur at a rate of unity, 
independent of the neighboring sites. The reverse transition, a vacant site becoming 
occupied, is only possible if at least one of its nearest neighbors (NNs) is occupied: 
the transition from crj = to <jj = 1 occurs at rate Ar, where r is the fraction of NNs 
of site i that are occupied. Thus the state cr^ = for all i is absorbing. At a certain 
critical value A c the system undergoes a phase transition between the active and 



the absorbing state [12] . The CP has been studied intensively via series expansion 
and Monte Carlo simulation, and its critical properties are known to high precision 

We now define a two-species symbiotic contact process (CP2S). Let the indicator 
variables for occupation of site i by species A and B be and rji, respectively. The 
allowed states for a site, (ai,r)i), are (0,0), (0,1), (1,0), and (1,1). The transitions 
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(0,0) — > (0,1) and (1,0) — > (1,1) occur at rate Xta with ta the fraction of NNs 
bearing a particle of species A. Similarly, the rate for the transitions (0, 0) — > (1, 0) 
and (0, 1) — > (1, 1) is Xr B with r B the fraction of NNs bearing a particle of species 
B. The transitions (0,1) — > (0,0) and (1,0) — > (0,0) occur at a rate of unity, 
whereas (1,1) —> (1,0) and (1,1) —> (0,1) at rate /i. This set of transition rates 
describes a pair of contact processes inhabiting the same lattice. For fi = 1 the two 
processes evolve independently, but for \i < 1 they interact symbiotically, since the 
annihilation rates are reduced at sites with both species present. We note that the 
rates are symmetric under exchange of species labels A and B. 

We also study a CP with creation at second-neighbor sites. In Ref. [l^] a modified 
CP was defined so: 

i) In addition to creation at NNs, at rate Ai, we allow creation at second neighbors, 
at rate A2. For bi-partite lattices such as the ring or the square lattice, Ai is the rate 
of creation in the opposite sublattice, while A2 is the rate in the same sublattice as 
the replicating particle. 

ii) The annihilation rate at a given site is 1 + un 2 , with n denoting the number 
of occupied NNs. 

For v > 0, the presence of particle in one sublattice tends to suppress their 
survival in the other, leading to the possibility of sublattice ordering, as discussed 
in Q. 

Suppose now that Ai = 0, and let A2 = A. Then the populations in the two sub- 
lattices constitute distinct species, since creation is always in the same sublattice. 
For v < 0, moreover, the two species interact in a symbiotic manner, analogous to 
that in the two-species CP defined above. (For v = the two sublattices evolve 
independently.) We call this process the symbiotic sublattice contact process (CP- 
SLS). 

Both the CP2S and CPSLS possess four phases: the fully active phase (nonzero 
populations of both species), a symmetric pair of partly active phases (only one 
species present), and the inactive phase (all sites inactive). The latter is absorbing 
while the partly active phases represent absorbing subspaces of the dynamics. (That 
is, a species cannot reappear once it goes extinct.) Let A Ci0 denote the critical 
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creation rate of the basic CP. In the CP2S with u = 1 (or the CPSLS with v = 0), 
the critical creation rate must be A c0 . The same applies for the transitions from 
the partly active phases to the absorbing one, regardless of the value of a or v. 
Intuitively, in the presence of symbiotic interactions, one expects the transition 
from the fully active to the absorbing phase to occur at some A c < A Cj o, since the 
annihilation rate is reduced. Since this expectation is borne out numerically, the 
partly active phases are of little interest, as they are not viable in the vicinity of the 
fully active-absorbing phase transition. Understanding the latter transition is the 
principal objective of this study. 

As a first step in characterizing the phase diagrams of the models, we develop 
mean-field approaches. The derivation of a dynamic mean-field theory (MFT) for 
an interacting particle system begins with the master equation for the set of one- 
site probabilities (or, more generally, the n-site joint probability distribution) 
In this equation, the n-site probability distribution is inevitably coupled to the 
distribution for n + 1 or more sites. An n-site MFT is obtained by estimating the 
latter distribution(s) in terms of that for n sites. Here we consider the simplest 
cases, n — 1 and 2. 

Consider the CP2S in the one-site approximation. Denoting the probabilities for 
a given site to be vacant, occupied by species A only, by species B only, and doubly 
occupied by p , pa, Pb, and Pab, respectively, assuming spatial homogeneity, and 
factorizing two-site joint probabilities (p[(cri, r]i) , (ffj,Vj)] = p[(°~«; Vi)}p[( a j, Vj)}) one 
readily obtains the equations 



dpo 
dt 
dpA 
dt 
dp B 
dt 

dpAB 



-Ap (PA + Pb) + Pa + Pb, 
ApoPa + Was - (1 + Ap B )pA, 
A^oPb + Wab ~ (1 + >^Pa)pb, 

KpaPb + PbPa) - 2fip A B, (1) 



dt 

where Pa = Pa + Pab and ps = Pb + Pab- If one species is absent (so that, say 
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Pb = Pab = 0) this system reduces to the MFT for the basic contact process, 
Pa = Ap^(l — Pa) — Pa, with a critical point at A = 1. To study the effect of 
symbiosis we seek a symmetric solution, Pa = Pb = P- in this case one readily finds 
the stationary solution: 



P- 

P 



2A(1 -fj) 
and 



2(1- A t)-A+ VA 2 -4/i(l-/i)l . (2) 



Ap 2 



For /i > 1/2, p grows continuously from zero at A = 1, marking the latter value as 
the critical point. The activity grows linearly, p ~ \p/(2p,— 1)](A — 1), in this regime. 
For fj, < 1/2, however, the expression is already positive for A = a/4/i(1 — yu) < 1, 
and there is a discontinuous transition at this point. The value fi — 1/2 may be 



viewed as a tricritical point; here p ~ VA — 1 for A > 1. Numerical integration of the 
MFT equations confirms the above results. For fi < 1/2, MFT in fact furnishes the 
spinodal values of A. For a given set of initial probabilities, the numerical integration 
converges to the active stationary solution for A > A* and to the absorbing state 
for smaller values of A. For the most favorable initial condition, i.e., Pab(0) — > 1, 
A* — > \(~^ = A/4yu(l — /i), the lower spinodal, while for a vanishing initial activity, 
Pa, Pb — > 0, A* — > = 1. The stationary activity at A* is nonzero. Fig. [1] shows 
the stationary probabilities versus A for fi = 1/4. 

The two-site MFT for the one-dimensional CP2S involves ten pair probabilities 
and a set of thirty-two transitions. The resulting phase diagram is qualitatively 
similar to that of the one-site MFT. For /i > 0.75, the transition is continuous and 
occurs at A = 2, the same value as for the basic CP at this level of approximation. 
There is a tricritical point at fi = 0.75, below which the transition is discontinuous; 
Fig. |2] shows the phase phase boundary. 

The one-site mean-field theory (MFT) for the CPSL was developed in Ref. fl^ |. 
Adapted to the present case (creation only in the same sublattice, symbiotic inter- 
action), the equation is 
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FIG. 1: Density p of species A (lower curve) and of doubly occupied sites, pab (upper curve) in 
the one-site approximation for the CP2S, /i = 0.25. 



^ = -(1 - vq 2 P 2 B ) PA + A 2PA (1 - p A ) (4) 

and similarly for pa ^ Pb, on a lattice of coordination number q. (Here pj denotes 
the fraction of occupied sites in sublattice j.) As we seek a symmetric solution, 
we set pa = Pb- The resulting equation yields a continuous phase transition at 
A = 1, independent of v. (Note that v must be greater than —1/16; smaller values 
correspond to a negative annihilation rate, for p near unity.) The two-site approxi- 
mation is likely to provide a better description of the CPSLS, since in this case the 
nearest-neighbor double occupancy probability is an independent variable, analo- 
gous to "Pab in the one-site MFT of the CP2S. Since such an analysis is unlikely to 
result in additional insights, we shall not pursue it here. 

Although MFT predicts a discontinuous phase transition in the CP2S in any 
number of dimensions, such a transition is not possible in one-dimensional systems 



with short-range interactions and free of boundary fields 



In one dimension the 



active-absorbing transition should be continuous, as we have indeed verified in sim- 



S 




FIG. 2: Phase boundary in the A — [i plane as given by two-site MFT for the CP2S on the line. 
The curved portion represents the lower spinodal, X^(n). 

ulations. Although our simulations show no evidence of a discontinuous transition 
in two dimensions (d = 2), such a transition remains a possibility for d > 2, for 
small values of /i. A discontinuous transition might also arise under rapid particle 
diffusion, as this generally favors mean-field- like behavior. 

III. SIMULATIONS 

We performed extensive Monte Carlo simulations of the CP2S on rings and on 
the square lattice (with periodic boundaries), and of the CPSLS on rings. A general 
observation is that both models appear to be more strongly affected by finite-size 
corrections than is the basic CP. 

In the simulation algorithm for the two-species CP, we maintain two lists, i.e., of 
singly and doubly occupied sites. Let N s and N d denote, respectively, the numbers 
of such sites, so that N p — N s + 2Nd is the number of particles. The total rate of 
(attempted) transitions is \N P + N s + 2fiNd = 1/At, where At is the time increment 
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associated with a given step in the simulation. At each such step, we choose among 
the events: (1) creation attempt by an isolated particle, with probability XN s At; (2) 
creation attempt by a particle at a doubly occupied site, with probability 2\NdAt; 
(3) annihilation of an isolated particle, with probability N s At; and (4) annihilation 
of a particle at a doubly occupied site, with probability 2/iiVd. Once the event type 
is selected we choose a site i from the appropriate list. In case of annihilation, a 
particle is simply removed, while creation requires the choice of a neighbor, j, of 
site i, and can only proceed if j is not already occupied by a particle of the species 
to be created. For creation by a particle at a doubly occupied site, the species of 
the daughter particle is chosen to be A or B with equal probability, and similarly 
for annihilation at a doubly occupied site. 

In simulations of the CPSLS we maintain a list of occupied sites. At each step 
a site is selected from the list; an attempt to create a new particle, at one of the 
second-neighbor sites, is chosen with probability p = A/(l + A2 + ^n\)\ the site is 
vacated with the complementary probability, 1 — p. The time increment associated 
with each event is At = 1/N P , with N p the number of particles just prior to the 
event. 



A. Results: CP2S in one dimension 

We studied the CP2S using three values of /1: 0.9, 0.75, and 0.25. While the 
first case may be seen as a relatively small perturbation of the usual CP (/x = 1), 
the third represents a very strong departure from the original model. We perform 



three kinds of studies: quasi-stationary (QS) 15] . initial decay (starting from a 
maximally active configuration), and spreading, in which the initial condition is a 
doubly occupied site in an otherwise empty lattice. Although the critical value, 
A c (/i), can be estimated using each method, spreading simulations proved the most 
effective in this regard. 

In the QS simulations, we study system sizes 800, 1600, 3200, 6400, and 12800, 
with each run lasting 10 7 time units; averages and uncertainties are calculated over 
10 - 80 runs. We use three well established criteria to estimate the critical value: 
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(1) power-law dependence of the order parameter on system size, p ~ L Pl v ±^ and 

(2) of the lifetime, r ~ L z , as well as (3) convergence of the moment ratio, m p (L), 



to a finite limit, m c , as L — > oo 16J. Here m p = (p 2 )/(p) 2 - The order parameter 
is defined as the density of individuals, i.e., p = (Na + Nb)/L. A related quantity 
of interest is the density q of doubly occupied sites; the moment ratio m q is defined 
in a manner analogous to m p . Two further quantities of interest are the scaled 
variances of p and q; we define Xp = L d vai(p) and similarly for Xq- The expected 
critical behavior is x ~ L 1 ^ 1 - , where the critical exponent 7 satisfies the hyperscaling 
relation 7 = du^ — 2/3 

A preliminary estimate of A c is obtained from the crossings of m p for successive 
system sizes, L and 2L. For p = 0.75, for example, this yields A c = 3.0337. The 
plot of m p and m q (see Fig. [3]) indicates that A c > 3.0336 (since m p curves upward), 
while the slight downward curvature for A = 3.0037 suggests that this value may 
be slightly above critical. This graph also suggests that m p and m q approach the 
same limiting value, despite marked differences for smaller system sizes. Table [I] 
summarizes our findings for the critical parameters obtained from QS simulations. 



TABLE I: CP2S in one dimension: results from QS simulations, L = 800, 1600, 3200, 
6400, and 12800. For p = 0.25 the maximum size is 6400. 





A c 




z 


m p 


m q 


(l/v±)p 




0.9 
0.75 
0.25 


3.2273(1) 
3.03370(5) 
1.76297(1) 


0.25(2) 
0.241(6) 
0.248(3) 


1.50(5) 
1.64(5) 
1.56(4) 


1.168(12) 
1.163(10) 
1.168(3) 


1.164(4) 
1.166(2) 
1.169(3) 


0.627(20) 
0.528(6) 
0.500(1) 


0.474(7) 
0.486(1) 
0.492(2) 


CP/DP 


3.29785 


0.25208(5) 


1.5807(1) 


1.1736(1) 




0.49584(9) 





The initial decay studies use, as noted above, an initial configuration with all sites 
doubly occupied. The activity then decays, following a power law, p ~ t~ 5 , at the 
critical point jl^], until it saturates at its QS value. The larger the system size, the 
longer the period of power-law decay, and the more precise the resulting estimate 
for the critical exponent 6; here we use L = 25600 or 51200. Averages are calculated 
over 500-3000 realizations. As the order parameter decays, its fluctuations build 
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FIG. 3: (Color online) QS Simulation of one-dimensional CP2S: moment ratios m p (filled symbols) 
and TTiq (open symbols) versus 1/L for the one-dimensional model with fj, = 0.75. The upper curve 
in each pair is for A = 3.0336, the lower for A = 3.0337. 

up; at the critical point, the moment ratio is expected to follow m — 1 ~ t l l z [isj]. 
Since we expect p and q to scale in the same manner, we define exponents 5 p and 
5 q , and, similarly, z p and z q , based on the behavior of m p and m q , respectively. 
Figure EB for p = 0.9, shows that p and q decay in an analogous manner, and 
follow power laws at long times, although there are significant deviations from a 
simple power law at short times; the decay exponents are consistent with the value 
of S for directed percolation in one space and one time dimension (see Table HT1) . 
The growth of fluctuations follows a more complicated pattern, as shown in Fig. 
[5j At relatively short times, m p — 1 ~ t 1 ^ Zp , with z p = 1.63(2), not very different 
from the DP value; m q — 1 also grows as a power law in this regime, but with an 
apparent exponent of z q = 2.06(1). At longer times z p appears to take a smaller 
value (1.31(1) for 7.5 < hit < 10.5), while z q shifts to a value close to that of DP 
(1.61(1) for 10 < hat < 14). The reason for the distinct behaviors of m p and m q , in 
marked contrast with the similar scaling of p(t) and q(t), is unclear. While scaling 
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anomalies are observed in the initial decay studies for fi = 0.9 and 0.75, for strong 
symbiosis (/i = 0.25) they are absent, as seen in Table Ull which summarizes the 
results of the initial decay studies. (In this table, the values listed for z p and z q 
reflect the latter part of the evolution, during which the order parameter decays in 
the expected manner.) 




-1 

c 

-2 

cr 

-3 
-4 

2 4 6 8 10 12 14 16 
Int 



FIG. 4: (Color online) Initial-decay simulation of CP2S in one dimension: decay of the particle 
density p (upper curve) and the density q of doubly-occupied sites in initial decay studies with 
/i = 0.9, A = 3.2273, system size L = 51200. The slopes of the regression lines are -0.161 (p) and 
-0.162 (q). 



TABLE II: CP2S in one dimension: results from initial-decay studies. 



n 


L 


A c 


5 


z p 




0.9 


51200 


3.2273 


0.161(1) 


1.44(1) 


1.54(2) 


0.75 


25600 


3.0337 


0.1625(10) 


1.48(4) 


1.55(4) 


0.25 


51200 


1.76297 


0.1581(3) 


1.56(1) 


1.58(1) 


DP 






0.1599 


1.5807(1) 
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FIG. 5: (Color online) Initial-decay simulation of CP2S in one dimension: growth of fluctuations 
in p (lower curve) and q (upper curve) for parameters as in Fig. 2) The slopes of the regression 
lines are (left to right): 0.613, 0.694, and 0.649. 

In the spreading studies, each realization runs to a maximum time of t m (unless 
it falls into the absorbing state prior to this). The system size is taken large enough 
so that activity never reaches the boundary. Here we use t m = 2 x 10 6 and L = 10 5 ; 
averages are calculated over 10 4 or 2 x 10 4 realizations. At the critical point, one 
expects to observe power-law behavior of the survival probability, P(t) ~ t~ s , the 
mean number of particles, n(t) ~ fj, and the mean-square distance of particles 
from the initial seed, R 2 (t) ~ t Zs 17|. Here 5 is the same exponent as governs 
the initial decay of the activity, and z s is related to the dynamic exponent z via 
z s = 2/z. Deviations from asymptotic power laws, indicating off-critical values of 
the control parameter A, are readily identified in spreading simulations, leading to 
precise estimates for A c . 

The spreading behavior is characterized by clean power laws, as illustrated in 
Fig. O As this plot makes clear, the mean particle number, n p , and the mean 
number of doubly occupied sites, ri2, grow with the same critical exponent. Precise 
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estimates of the spreading exponents are obtained via analysis of local slopes such as 
S(t), defined as the inclination of a least-square linear fit to the data (on logarithmic 
scales), on the interval [t/a, at}. (The choice of the factor a represents a compromise 
between high resolution, for smaller a, and insensitivity to fluctuations, for larger 
values; here we use a = 4.59.) Curvature in a plot of a local slope versus 1/t signals 
an off-critical value. Figure [7] shows the behavior of S(t) for \i = 0.25. The spreading 
exponents, summarized in Table IHH are in good agreement with the values for DP 
in 1+1 dimensions. (We note that in all three cases, r] p = rj 2 to within uncertainty.) 




FIG. 6: (Color online) Spreading simulation of CP2S in one dimension: survival probability P(t), 
total particle number n p (t), number of doubly occupied sites, Ti2(t), and mean-square distance 
from seed, R 2 (t). Parameters fx = 0.25, A = 1.76297. 



B. Contact process with creation at second-neighbors 

We studied the CPSLS using QS and initial decay simulations. The results from 
the former, based on FSS analysis of studies using L = 800, 1600, 3200, 6400, 
and 12800, are summarized in Table IIVI The value of uj_ was estimated (for v = 

15 



-0.15 




-0.19 ' 1 1 1 1 ^ ^ 

0.00000 0.00001 0.00002 0.00003 

1/t 



FIG. 7: (Color online) Spreading simulation of CP2S in one dimension: local slope S(t) versus 1/t 
for \i = 0.25 and (lower to upper) A = 1.7629, 1.76295, 1.76297, and 1.7630. 

TABLE III: One-dimensional CP2S: Results from spreading simulations. 





A c 


5 


V 


Z s 


0.9 


3.2273 


0.165(1) 


0.310(1) 


1.257(2) 


0.75 


3.0337 


0.1595(5) 


0.3180(5) 


1.265(5) 


0.25 


1.76297 


0.158(1) 


0.315(3) 


1.265(10) 


DP 


3.29785 


0.15947(5) 


0.31368(4) 


1.26523(3) 



—0.1 only) via analysis of the derivatives \dm/d\\, dlnr/dX and dlnp p /dX in the 
neighborhood of the critical point. Finite-size scaling implies that the derivatives 
follow \dx/dp\ oc L 1 ^ 1 - (here x stands for any of the quantities mentioned). We 
estimate the derivatives via least-squares linear fits to the data on an interval that 
includes A c . (The intervals are small enough that the graphs show no significant 
curvature.) Linear fits to the data for m, lnp p , and lnr yield 1/Vj_ = 0.94(2), 
0.92(3), and (again) 0.92(3), respectively, leading to the estimate u± = 1.08(3). 
Results of the initial decay studies are summarized in Table [V] As in the two- 
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species CP, the value of z obtained from analysis of m(t) appears to be smaller than 
the DP value, whereas the result obtained from QS simulations is consistent with 
that of DP. 

TABLE IV: One-dimensional CPSLS: results from quasistationary simulations. 



V 


A c 




z 


m c 


V± 


-0.05 
-0.1 
-0.2 


3.1489(1) 
2.8878(1) 
2.0502(1) 


0.235(8) 
0.242(1) 
0.253(6) 


1.63(5) 
1.612(12) 
1.59(1) 


1.154(5) 
1.161(3) 
1.170(6) 


1.08(3) 


DP 


3.29785 


0.25208(5) 


1.5807(1) 


1.1736(1) 


1.096854(4) 



TABLE V: One-dimensional CPSLS: results from initial decay simulations. 



V- 


L 


A c 


5 


z 


-0.05 


50000 


3.1489 


0.1458(5) 


1.45(2) 


-0.1 


50000 


2.8878 


0.1484(7) 


1.45(3) 


-0.2 


20000 


2.0503 


0.1597(3) 


1.53(1) 


DP 






0.1599 


1.5807(1) 



C. Two-species contact process in two dimensions 

We performed extensive Monte Carlo simulations of the CP2S on square lattices 
using both initial decay and quasi- stationary (QS) simulations. In order to locate the 
critical point with good precision, we study the initial decay of the particle density, 
starting from a maximally active initial condition (all sites doubly occupied). We 
use lattices of linear size L = 4000, and average over at least 20 different realizations. 
Figure |8] shows the decay of p(t) for \x = 0.1. After an initial transient, during which 
the density evolves slowly, the particle density follows a power law with 5 = 0.46(1), 
compatible with the value (5 = 0.4523(10)) for the DP class in 2+1 dimensions. 
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The transient behavior lasts longer, the larger is /x, as shown in Fig. [9j However the 
relaxation is seen to cross over to DP-like behavior for all values studied, except for 
fi = 0.9, for which the transient regime persists throughout the entire simulation. 
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FIG. 8: (Color online) CP2S in two dimensions: Density of active sites starting from a maximally 
active initial condition, for fi = 0.1, and A values as indicated. System size: L = 4000. 
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FIG. 9: (Color online) CP2S in two dimensions: density of active sites starting from a maximally 
active initial condition, for /i = 0.1, 0.25, 0.5, 0.75 and 0.9 (from top to bottom) and A = A c (/x) (see 
Table ED- The slope of the dashed line is -0.45. System size: L = 4000. 

Having determined A c to good precision in the initial decay studies, we perform 
QS simulations of the model on square lattices of linear size L = 20, 40, 320 with 
periodic boundaries. Figure [10] shows moment-ratio crossings and the finite-size 
scaling behavior of the density and lifetime for /i = 0.1. For the larger sizes we 
obtain (3/v± = 0.78(1) and z = 1.74(2), in good agreement with the best estimates 
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TABLE VI: Simulation: critical parameters for the two-dimensional CP2S. 



n 


A c 


P/v± 


z 


5 


m p 


m q 


0.9 


1.64515(5) 


0.63(5) 


1.95(5) 


> 0.35 


1.40(2) 


1.52(2) 


0.75 


1.61640(5) 


0.73(5) 


1.78(6) 


0.44(3) 


1.32(3) 


1.33(3) 


0.5 


1.47290(5) 


0.74(3) 


1.72(3) 


0.46(2) 


1.298(8) 


1.322(8) 


0.25 


1.13730(5) 


0.76(2) 


1.73(2) 


0.45(2) 


1.30(2) 


1.31(2) 


0.1 


0.743160(5) 


0.78(1) 


1.73(2) 


0.46(1) 


1.305(10) 


1.315(12) 


CP/DP 


1.64874(4) 


0.797(3) 


1.7674(6) 


0.4523(10) 


1.3264(5) 





for DP in 2+1 dimensions. Simulation results for the two-dimensional model are 
summarized in Table I VI I 

2.0 
1.8 
1.6 
1.4 
1.2 

'o°740 0.742 0.744 0.746 0.748 0.750 

X 

FIG. 10: (Color online) CP2S in two dimensions: QS moment ratio of particles vs. A. , for fj, = 0.1 
(System sizes: L = 40, 80, 160, 320 in order of steepness). Inset: QS density of active sites (circles), 
density of doubly occupied sites (squares) and lifetime of the QS state (triangles), for fj, = 0.1. 




IV. FIELD THEORETIC ANALYSIS 



In this section we extend the field theory or continuum representation of DP to 
the two-species case, to determine whether the presence of aditional species changes 
the scaling behavior. Since the theory of DP has been known for some time, we 
give a bare outline of this analysis, referring the reader to references 19N23I] for 
details. To begin, we modify the lattice model so as to facilitate the definition of 
a continuum description following the Doi-Peliti formalism 
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applied to DP in 



20|, which has been 



2l| and 23j . (The latter study applies the Wilson renormalization 
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group to the problem.) 

In the Doi-Peliti formalism, the master equation governing the evolution of the 
probability vector \P(t)) = (the sum is over all configurations), is 

written in the form d\P)/dt = L\P), where the evolution operator L is composed of 
creation and annihilation operators. Starting from this "microscopic" description, 
one derives an effective action S via a path-integral mapping. Then, taking the 
continuum limit, one arrives at a field theory for the model. Of the many lattice 
models that belong to the DP universality class, the simplest to analyze in this 
manner is the Malthus-Verhulst process (MVP). Here, each site i of a lattice hosts a 
number rii > of particles. The transitions at a given site are creation (rii — > n± + 1) 
at rate Xrii and annihilation [rii — > ri\ — 1) at rate \LUi + vn^rii — 1). In addition, 
particles hop between nearest-neighbor sites at rate D. 

For the MVP on a ring of I sites, one has the set of basis configurations \ni, ...,rie). 
Letting Cj and c] denote, respectively, annihilation and creation operators associated 
with site i, we have, by definition, Cj|ni, rii, ng) = n,|ni, n, — l,...,n^) and 
c] | rii, rii, ...,ne) = \rii, rii + 1, ng). Then the evolution operator for the MVP 
is. 



Lmvp = ^ [a(c] - 1)c|q + (1 - c\)(fi + vci)c { 

i 



(5) 



Following the steps detailed in 



21| . one arrives at the effective action for the MVP, 



Smvp = dt dx 



ip(d t + w- DV' z )^p + v^ 2 - Xftip 



(6) 



where w = fi — A, the continuum limit has been taken, and terms higher than third 
order have been discarded, as they are irrelevant to critical behavior. (We recall that 
i[)(x, t) is an auxiliary field that arises in the mapping. The operator that governs the 
evolution of the probability generating function is given by the functional integral 
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U t = JVipJ Xtyexp[-S(V>,$)]; see |2o|, 121 1 . 



The action of Eq. ([6]) is equivalent that of DP, and serves as the starting point for 
renormalization group (RG) analyses 0, lo[ 2s|. (One usually imposes the relation 
v = X via a rescaling of the fields, but this is not needed here.) In the RG analysis the 
bilinear term, naturally, defines the propagator, while the cubic terms correspond 
to the vertices shown in Fig. [TTJ These terms lead, via diagrammatic analysis, to 
a nontrivial DP fixed point below d c — 4 dimensions. The one-loop diagrams which 
yield, to lowest order, the recursion relations for parameters w, A, and v are shown 
in Fig. [U 





FIG. 11: The two three-field vertices in the field theory of DP. Lines exiting a vertex correspond 
to ip while those entering correspond to ip. 



Now consider the two-species CP. To formulate a minimal field theory, we consider 
a two-species MVP; call it MVP2. Let rrtj and n.j denote, respectively, the number 
of particles of species A and B at site i, and let aj and a\, and 6j and b\, denote 
the associated annihilation and creation operators. We require the annihilation rate 
for species A to be a decreasing function of and vice-versa; a simple choice for 
the annihilation rate of an A particle at site % is /iexp[— 77^], where 7 is a positive 
constant, and similarly for B particles, with rij replaced by r?v This corresponds to 
the evolution operator: 



21 



FIG. 12: The one-loop diagrams in the field theory of DP, leading to renormalization of w, fi, and 
A, respectively. 



'MVP1 



tV,, -76itfei 



D 



^ [A(6j - !)&}&, + (1 - m^ a ^ + vtybi 



(7) 



To avoid ambiguity, we interpret the exponentials as being in normal order, i.e., all 
creation operators to the left of annihilation operators. Recalling that terms with 
four or more fields are irrelevant, we may expand the exponentials, retaining only 
the terms oc b\bi and oc a\ai. Using : X : to denote the normal-ordered expression 
of X, it is straightforward to show that 



: e 



-7&t& ._ 



1 - (1 -e-^tfb + I = 1 



(8) 



where I consists of terms with four or more operators. (With the truncation comes 
the possibility of a negative rate, but this is of no consequence in the RG analysis.) 
Now, following the usual procedure, we obtain the effective action for the two-species 
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MVP: 




= j dt j dx ip{d t + w- DV 2 }il)) + vi)i? - h\) 2 ^) 
+ J dt J dx [0(d t +w- £>V 2 )y?) + vipip 2 - \0 2 tp] 
— v dt dx (pipip + ijjipLp , 



(9) 



where v — jfj,. Here ip and ip are fields associated with species A; ip and y3 are 
associated with species B. The first two lines of the above expression correspond to 
independent MVPs; the third represents the symbiotic interaction between them. 
[While such a minimal action could have been "postulated" directly, we prefer to 
start with the microscopic expression of Eq. ([7]), since it describes a valid stochastic 
process.] 

There are two cubic terms in the action involving only species A (i.e., the vertices 
shown in Fig. ITT|) . two involving only B (those of Fig. [11] drawn, say, with broken 
lines) and two vertices with a mixed pairoi incoming lines, and a single outgoing line, 
which may belong to either species. One readily identifies the one-loop diagrams 
leading to renormalization of the parameter v. On the other hand, no diagrams (at 
any order) involving mixed-species vertices can affect the recursion relations for the 
DP parameters w, u, and A. The reason is that the presence of a mixed-species 
vertex anywhere in a diagram implies that the lines entering the diagram are mixed, 
so that it can only contribute to the recursion relation for v. We conclude that the 
interaction between species cannot alter the scaling behavior, which must therefore 
remain that of DP. At one-loop order, there are two fixed-point values for u, namely, 
2A and zero, the latter corresponding to independent processes. 

V. CONCLUSIONS 

We study symbiotic interactions in contact-process-like models in one and two 
dimensions. For this purpose, we propose a two-species model (CP2S), in which 
the death rate is reduced (from unity to fi), on sites occupied by both species. A 
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related model (CPSLS), in which each species is confined to its own sublattice, is 
also studied in one dimension, and found to exhibit similar behavior. Simulations 
reveal that the phase transition between active and absorbing states is continuous, 
and that the critical creation rate A c is reduced in the presence of symbiosis. This 
means that the loss of one species will rapidly lead to extinction, since the system 
is then a basic contact process operating at A < A c . Although this might suggest 
identifying the density q of doubly occupied sites as the order parameter, we find 
that the particle density p (which includes a large contribution from singly occupied 
sites), scales in the same manner as q. 

Mean-field theory (in both the one- and two-site approximations), predicts a 
discontinuous phase transition in any number of dimensions, for p sufficiently small. 
A discontinuous transition between an active and an absorbing phase is not expected 
in one-dimensional systems of the kind studied here 24J, nor do our simulations 
show any evidence of a discontinuous transition in two dimensions. Nevertheless, 
we cannot discard the possibility of such a transition for d > 2, for small values of 
p, or under rapid particle diffusion, which generally favors mean-field-like behavior. 

Overall, the critical behavior of the symbiotic models is consistent with that of 
directed percolation. Corrections to scaling are, however, more significant than in 
the basic CP, so that a study restricted to smaller systems, or to only one kind of 
simulation, could easily suggest non-DP behavior. These corrections are stronger, 
and of longer duration, the smaller the intensity of symbiosis. Thus, in the two- 
dimensional case, the decay of p (in initial-decay studies) attains the expected power- 
law regime (with a DP value for the decay exponent), except for p = 0.9, the weakest 
symbiosis studied. A similar tendency is observed in the QS simulations of the one- 
dimensional CP2S, for which the estimates for critical exponents and the critical 
moment ratio m c differ most from DP values for p = 0.9. 

In the initial-decay studies in one dimension, for smaller intensities of symbiosis 
(i.e., p = 0.9 and 0.75), we observe anomalous growth of fluctuations in the order 
parameter. The latter are characterized by m p — 1 = var(p)/p 2 , which is expected to 
grow ~ t 1//z , before saturating at its QS value. The growth at long times corresponds 
to a z value significantly smaller than that of DP. The exponent z q associated with 



24 



the growth of m q is substantially larger, though still slightly below the DP value. In 
contrast with these anomalies, the spreading exponents are found to take DP values 
in one dimension, independent of the degree of symbiosis. Thus we are inclined to 
regard the asymptotic scaling of the symbiotic models as being that of DP, and to 
interpret the deviations as arising from finite-time and finite-size corrections. One 
might conjecture that under strong symbiosis, the critical system is rapidly attracted 
to the DP fixed point (although not as rapidly as is the basic CP), whereas for weak 
symbiosis, it makes a long excursion into a regime in which DP-like scaling is not 
evident, before finally returning to the vicinity of the DP fixed point. The asymptotic 
scaling behavior is presumably associated with large, sparsely populated but highly 
correlated regions of doubly occupied sites, which, for reasons of symmetry, behave 
analogously to DP space-time clusters. The presence of isolated particles, which are 
relatively numerous and long-lived for weak symbiosis, could mask the asymptotic 
critical behavior, on short scales. We defer further analysis of these questions to 
future work. 

Extending the field theory of DP to the two-species case, we find that the ir- 
relevance of four- field terms makes DP extremely robust, since the only possible 
three-field vertices are already present in the single-species theory. This means that 
the interaction between species cannot alter the scaling behavior, as already noted 



by Janssen in the case of multi-species DP processes [13|] . Our simulation results, as 
noted, support this conclusion. A more detailed field-theoretic analysis, including 
the evolution of the lowest order irrelevant terms, might shed some light on the 
scaling anomalies observed in the simulations. 
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